Esta guía explica cómo realizar algunas tareas comunes de análisis de datos de recuento de secuenciación microbiana utilizando tidytacos. Ilustraremos utilizando un conjunto de datos con muestras de microbioma humano del tracto vaginal, tomadas de este artículo por Lebeer et al.
Un objeto tidytacos es simplemente una lista de tres tablas:
counts: estos son los recuentos de lecturas para cada taxón (OTU/ASV/filotipo) en cada muestra. Cada fila representa dicho recuento de lecturas.
samples: esta tabla contiene los metadatos de muestra. Cada fila representa una muestra.
taxa: esta tabla contiene la taxonomía y otros metadatos de los taxones. Cada fila representa un taxón.
El paquete se llama tidytacos porque cada una de las tablas está ordenada (en inglés tidy): cada fila representa una observación y cada columna una variable (puedes encontrar más información sobre la ordenación de datos en este enlace).
En caso de que aún no hayas instalado tidytacos, puedes instalarlo usando devtools:
install.packages("devtools")
devtools::install_github("LebeerLab/tidytacos")
Para esta guía, sólo necesitamos cargar tres paquetes: tidytacos (por supuesto), el conjunto de paquetes tidyverse, y VTutorials.
library(tidyverse)
library(tidytacos)
library(VTutorials)
Nuestro conjunto de datos de ejemplo está disponible en el paquete VTutorials y no es necesario importarlo ni convertirlo. Se llama “vdata” y comenzamos inspeccionando la tabla de muestras:
glimpse(vdata$samples)
## Rows: 200
## Columns: 28
## $ sample <chr> "SAMEA13411869", "SAMEA13410951",…
## $ sample_id <chr> "s16", "s100", "s114", "s170", "s…
## $ Technical.run_20201009_isala_cross_1 <lgl> FALSE, TRUE, FALSE, FALSE, TRUE, …
## $ Technical.run_20201016_isala_cross_2 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201023_isala_cross_3 <lgl> TRUE, FALSE, TRUE, FALSE, FALSE, …
## $ Technical.run_20201030_isala_cross_4 <lgl> FALSE, FALSE, FALSE, TRUE, FALSE,…
## $ Technical.run_20201106_isala_cross_5 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201120_isala_cross_6 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201204_isala_cross_8 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201211_isala_cross_9 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20210119_isala_cross_7 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.libsize <dbl> 37957, 22782, 32150, 7707, 24123,…
## $ General.Age <dbl> 24, 60, 24, 33, 29, 24, 23, 42, 3…
## $ Health.BMI <dbl> 30.11938, 19.00391, 25.03414, 23.…
## $ Health.Antibiotic.3months <lgl> FALSE, TRUE, TRUE, TRUE, FALSE, F…
## $ Sexual.Intercourse.24hours <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ plate <chr> "PlateBA", NA, "PlateBB", "PlateB…
## $ volume <dbl> 2.2025248, NA, 1.5580467, 1.83601…
## $ dna_conc <dbl> 11.169, NA, 15.789, 20.648, 24.65…
## $ read_conc <dbl> 17233.4038, NA, 20634.8110, 4197.…
## $ read_conc_norm <dbl> 0.9693078, NA, 1.1606229, 0.57860…
## $ lib_size <dbl> 37957, NA, 32150, 7707, 24123, 29…
## $ `Sample month` <dbl> 7, NA, 7, 7, 7, 7, 7, 7, 7, 7, 7,…
## $ `Transit time` <dbl> 1, NA, 1, 1, 1, 2, 1, 1, 6, 2, 1,…
## $ valencia_subcst <chr> "I-A", "IV-C0", "III-A", "IV-B", …
## $ valencia_cst <chr> "I", "IV-C", "III", "IV-B", "IV-B…
## $ max_genus1 <chr> "L. crispatus", "Prevotella", "L.…
## $ max_genus2 <chr> "Gardnerella", "Other", "Prevotel…
Luego echamos un vistazo rápido al número total de muestras, ASV y lecturas (reads en inglés) en el objeto tidytacos:
tacosum(vdata)
## n_samples n_taxa n_reads
## 200 942 5022096
Podemos crear muy fácilmente un gráfico para explorar un subconjunto de nuestras muestras (por ejemplo, solo muestras de participantes que han tomado antibióticos en los últimos 3 meses y son menores de 40 años) de la siguiente manera:
vdata %>%
filter_samples(Health.Antibiotic.3months == TRUE, General.Age <= 40) %>%
tacoplot_stack()
La función filter_samples hace lo que dice: filtrar muestras. También
eliminará los taxones de la tabla de taxones que tengan cero lecturas
totales en las muestras restantes. La función tacoplot_stack devuelve
una buena visualización de diagramas de barras apiladas de los taxones
más abundantes en nuestras muestras.
¿Se te ocurren otras ideas de como filtrar y rápidamente visualizar vdata?
Nuestra siguiente pregunta para este conjunto de datos es hasta qué punto la actividad sexual dentro de las 24 horas antes de tomada las muestras afecta la composición microbiana de las participantes en edad reproductiva (asumamos <=40 años). Para hacernos una idea primero lo podemos visualizar:
vdata_less40 <- vdata %>% filter_samples(General.Age <= 40)
tacoplot_stack(vdata_less40)+
geom_point(aes(y=-0.02,color=Sexual.Intercourse.24hours))
## Diversidad Alpha Para explorar la diversidad alfa, creemos una
versión enrarecida (rarefied en inglés) del conjunto de datos:
vdata_rar <- vdata %>%
add_total_count() %>%
filter_samples(total_count >= 2000) %>%
rarefy(2000) %>%
add_alpha()
La función add_total_count agregará números totales de lectura de muestra a la tabla de muestra.
La función rarefy submuestreará aleatoriamente todas las muestras n veces. Solo funciona si el recuento de lecturas de cada muestra es igual o superior a n.
Para determinar la riqueza de ASV, optamos por enrarecer primero, pero esto puede depender de sus datos.
La función add_alpha se puede utilizar para agregar varias métricas de diversidad alfa a la tabla de muestra.
Podemos visualizar la diversidad alpha de muestras de participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no:
vdata_rar %>%
samples() %>%
ggplot(aes(x = Sexual.Intercourse.24hours, y = observed, fill = Sexual.Intercourse.24hours)) +
geom_boxplot(outlier.shape = NA)+
geom_jitter(height = NULL)
Nos gustaría abordar las diferencias entre muestras de participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no. También estamos más interesados en los géneros que en los ASV.
Una PCoA podría ofrecer información:
vdata_genus <- vdata %>%
aggregate_taxa(rank = "genus")
tacoplot_ord_ly(vdata_genus, Sexual.Intercourse.24hours, samplenames = sample,
dim = 3)
La función aggregate_taxa fusiona todas las filas de la tabla de taxones en un nivel taxonómico específico, en este caso el nivel de género. Como ocurre con todas las funciones de tidytacos, todas las demás tablas del objeto tidytacos se ajustan en consecuencia.
La función tacoplot_ord_ly determinará la abundancia relativa de taxones en las muestras y luego utilizará las diferencias de Bray-Curtis para ordenar muestras en un espacio bidimensional (o tridimensional) según su composición taxonómica. La adición argumental de plotly “_ly” hace que la figura sea interactiva, lo cual es realmente bueno para el trabajo exploratorio. Esto también funciona para otras funciones de visualización.
La siguiente pregunta lógica es hasta qué punto la actividad sexual reciente (dentro de las 24 horas) determina la variabilidad de la composición de la comunidad microbiana. No olvidemos que la composición microbiana vaginal varia con la edad, por eso incluyamos edad en el modelo.
perform_adonis(vdata_genus, c("General.Age", "Sexual.Intercourse.24hours"))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = as.formula(paste("counts_matrix", formula_RHS, sep = " ~ ")), data = metadata, permutations = permutations)
## Df SumOfSqs R2 F Pr(>F)
## General.Age 1 1.717 0.02676 5.5016 0.001 ***
## Sexual.Intercourse.24hours 1 0.956 0.01489 3.0607 0.018 *
## Residual 197 61.500 0.95835
## Total 199 64.173 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
La función perform_adonis realizará un ANOVA
PERMutacional para determinar el efecto de las variables de muestra en
las diferencias de Bray-Curtis de las comunidades. El resultado muestra
que la edad es un contribuyente a la composición de la comunidad
microbiana (R cuadrado = 0.02676).
A continuación, nos gustaría saber cuáles de los 20 géneros más abundantes son significativamente más abundantes en las participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no.
vdata_genus <- vdata_genus %>% add_codifab(Sexual.Intercourse.24hours,
max_taxa = 20)
vdata_genus$taxon_pairs <- filter(vdata_genus$taxon_pairs, wilcox_p < 0.05)
tacoplot_codifab(vdata_genus, FALSE_vs_TRUE)
La función
add_codifab agregará una tabla llamada
taxon_pairs al objeto tidytacos, con la abundancia diferencial del taxón
entre las dos condiciones (con respecto al taxón de referencia), para
cada par de un taxón y un taxón de referencia.
La función tacoplot_codifab devuelve un gráfico para
visualizar la abundancia diferencial de taxones entre condiciones, en
comparación con todos los demás taxones como referencia. Podemos
observar que es más probable que Staphiloccus y L. iners group sea
típico encontrar cuando una participante ha tenido actividad sexual
dentro de las 24 horas antes de tomada las muestras.
Es de destacar que existen muchos métodos de análisis de abundancia diferencial y ninguno de ellos es perfecto. Interprete tus resultados con cuidado.
Descargo de responsabilidad: Esta guía es una adaptación de la versión en inglés del tutorial de tidytacos.